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Abstract — We propose an image processing scheme based on 
reordering of its patches. For a given corrupted image, we 
extract all patches with overlaps, refer to these as coordinates 
in high-dimensional space, and order them such that they are 
chained in the "shortest possible path", essentially solving the 
traveling salesman problem. The obtained ordering applied to 
the corrupted image, implies a permutation of the image pixels 
to what should be a regular signal. This enables us to obtain 
good recovery of the clean image by applying relatively simple 
one-dimensional (ID) smoothing operations (such as filtering or 
interpolation) to the reordered set of pixels. We explore the use 
of the proposed approach to image denoising and inpainting, and 
show promising results in both cases. 

Index Terms — patch-based processing, traveling salesman, 
pixel permutation, denoising, inpainting. 



I. Introduction 

In recent years, image processing using local patches has 
become very popular and was shown to be highly effective - 
see [1-11] for representative work. The core idea behind these 
and many other contributions is the same: given the image 
to be processed, extract all possible patches with overlaps; 
these patches are typically very small compared to the original 
image size (a typical patch size would be 8 x 8 pixels). 
The processing itself proceeds by operating on these patches 
and exploiting interrelations between them. The manipulated 
patches (or sometimes only their center pixels) are then put 
back into the image canvas to form the resulting image. 

There are various ways in which the relations between 
patches can be taken into account: weighted averaging of 
pixels with similar surrounding patches, as the NL-Means 
algorithm does |1|, clustering the patches into disjoint sets 
and treating each set differently, as performed in Q, |[3|, ||4j 
|[5|, ij^, seeking a representative dictionary for the patches 
and using it to sparsely represent them, as practiced in (Tj, 
jSI, |9| and fTOl, gathering groups of similar patches and 
applying a sparsifying transform on them |TT| . A common 
theme to many of these methods is the expectation that every 
patch taken from the image may find similar ones extracted 
elsewhere in the image. Put more broadly, the image patches 
are believed to exhibit a highly- structured geometrical form in 
the embedding space they reside in. A joint treatment of these 
patches supports the reconstruction process by introducing a 
non-local force, thus enabling better recovery. 
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In our previous work p2| and p3| we proposed yet another 
patch-based image processing approach. We constructed an 
image-adaptive wavelet transform which is tailored to sparsely 
represent the given image. We used a plain ID wavelet 
transform and adapted it to the image by operating on a 
permuted order of the image pixelq^ The permutation we 
proposed is drawn from a shortest path ordering of the image 
patches. This way, the patches are leveraged to form a multi- 
scale sparsifying global transform for the image in question. 

In this paper we embark from our earlier work as reported in 
1 12 1 and |13|, adopting the core idea of ordering the patches. 
However, we discard the globality of the obtained transform 
and the processing that follows, the multi-scale treatment, and 
the sparsity-driven processing that follows. Thus, we propose 
a very simple image processing scheme that relies solely on 
patch reordering. We start by extracting all the patches of size 
^/n X ^/n with maximal overlaps. Once these patches are 
extracted, we disregard their spatial relationships altogether, 
and seek a new way for organizing them. We propose to refer 
to these patches as a cloud of vectors/points in R^, and we 
order them such that they are chained in the "shortest possible 
path", essentially solving the traveling salesman problem 
This reordering is the one we have used in |12| and 



but as opposed to our past work, our treatment from this 
point varies substantially. A key assumption in this work is 
that proximity between two image patches implies proximity 
between their center pixels. Therefore if the image mentioned 
above is of high-quality, the new ordering of the patches is 
expected to induce a highly regular (smooth or at least a 
piece- wise smooth) ID ordering of the image pixels, being 
the center of these patches. When the image is deteriorated 
(noisy, containing missing pixels, etc.), the above ordering is 
expected to be robust to the distortions, thereby suggesting 
a reordering of the corrupted pixels to "what should be" a 
regular signal. Thus, applying relatively simple ID smoothing 
operations (such as filtering or interpolation) to the reordered 
set of pixels should enable good recovery of the clean image. 

This is the core process we propose in this paper - for a 
given corrupted image, we reorder its pixels, operate on the 
new ID signal using simplified algorithms, and reposition the 
resulting values to their original location. We show that the 
proposed method, applied with several randomly constructed 
orderings and combined with a proposed subimage averaging 
scheme, is able to lead to state-of-the-art results. We explore 
the use of the proposed image reconstruction scheme to image 
denoising, and show that it achieves results similar to the 

^Note that the idea of adapting a wavelet transform to the image by 
reordering its pixels appeared already in 1 14], but the scheme proposed there 
did not use image patches, and targeted image compression only. 
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ones obtained with the K-SVD algorithm |T|. We also explore 
the use of the proposed image processing scheme to image 
inpainting, and show that it leads to better results compared to 
the ones obtained with a simple interpolation scheme and the 
method proposed in 1 16] which employs sparse representation 
modeling via the redundant DCT dictionary. Finally, we draw 
some interesting ties between this scheme and BM3D rationale 

m 

The paper is organized as follows: In Section II we in- 
troduce the basic image processing scheme. In Section III 
we explain how the performance of the basic scheme can be 
improved using a subimage averaging scheme, and describe 
the connection between the improved scheme and the BM3D 
algorithm. In Section IV we explore the use of the proposed 
approach to image denoising and inpainting, and present 
experimental results that demonstrate the advantages of the 
proposed scheme. We summarize the paper in Section IV with 
ideas for future work along the path presented here. 

II. Image Processing using Patch Ordering 

A. The Basic Scheme 

Let Y be an image of size A^i x N2 where A^iA^2 = 
and let Z be a corrupted version of Y, which may be noisy 
or contain missing pixels. Also, let z and y be the column 
stacked representations of Z and Y, respectively. Then we 
assume that the corrupted image satisfies 



permutations 



z My + V 



(1) 



where the N x N matrix M denotes a linear operator which 
corrupts the data, and v denotes an additive white Gaussian 
noise independent of y with zero mean and variance a^. In 
this work the matrix M is restricted to represent a point- 
wise operator, covering applications such as denoising and 
inpainting. The reason for this restriction is the fact that we 
will be permuting the pixels in the image, and thus spatial 
operations become far more complex to handle. 

Our goal is to reconstruct y from z, and for this end we 
employ a permutation matrix P of size N x N. We assume 
that when P is applied to the target signal y, it produces a 
smooth signal y^ = Py. We will explain how such a matrix 



may be obtained using the image patches in Section II-B We 
start by applying P to z and obtain z^ = Pz. Next, we take 
advantage of our prior knowledge that y^ should be smooth, 
and apply a "simple" ID smoothing operator H on z^, such 
as ID interpolation or filtering. Finally, we apply P~^ to the 
result, and obtain the reconstructed image 



(2) 



In order to better smooth the recovered image, we use an 
approach which resembles the cycle spinning method [17J . 
We randomly construct K different permutation matrices P/^, 
utilize each to denoise the image z using the scheme described 
above, and average the results. This can be expressed by 

1 ^ 

y = ^T.^k'H{p,z}. (3) 



k=l 



Fig. [T] shows the proposed image processing scheme. We next 
describe how we construct the reordering matrix P. 



ID inverse 
processing permutations 



xl/K 



Fig. 1: The basic image processing scheme. 



B. Building the Permutation Matrix P 

We wish to design a matrix P which produces a smooth 
signal when it is applied to the target image y. When the 
image Y is known, the optimal solution would be to reorder 
it as a vector, and then apply a simple sort operation on the 
obtained vector. However, we are interested in the case where 
we only have the corrupted image Z. Therefore, we seek a 
suboptimal ordering operation, using patches from this image. 

Let i/i and Zi denote the ith samples in the vectors y and z, 
respectively. We denote by the column stacked version of 
the ^/n X ^/n patch around the location of in Z. We assume 
that under a distance measure ( x^, Xj), proximity between 
the two patches x^ and Xj suggests proximity between the 
uncorrupted versions of their center pixels i/i and yj . Thus, we 
shall try to reorder the points x^ so that they form a smooth 
path, hoping that the corresponding reordered ID signal y^ 
will also become smooth. The "smoothness" of the reordered 
signal y^ can be measured using its total-variation measure 



Iy^IItf 



N 

E 

J=2 



(4) 



Let {x^}^^ denote the points {x^}^^ in their new order. 
Then by analogy, we measure the "smoothness" of the path 
through the points x^ by the measure 



N 



X^^ = ^Kx^,x^_,). 



(5) 



J=2 



Minimizing X^y comes down to finding the shortest path that 
passes through the set of points x^, visiting each point only 
once. This can be regarded as an instance of the traveling sales- 
man problem |15|, which can become very computationally 
expensive for large sets of points. We choose a simple approx- 
imate solution, which is to start from a random point and then 
continue from each point Xj^ to its nearest neighbor Xj^ with 



a probability pi ex exp 



, or to its second nearest 

w(XjqPCj2)' 



neighbor Xj2 with a probability p2 cx exp 
where e is a design parameter, and x^^ and Xj2 are taken 
from the set of unvisited points. 

We restrict the nearest neighbor search performed for each 
patch to a surrounding square neighborhood which contains 
B X B patches. When no unvisited patches remain in that 
neighborhood, we search for the nearest neighbor among all 
the unvisited patches in the image. This restriction decreases 

^Throughout this paper we will be using variants of the squared Euclidean 
distance. 
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(a) (b) 

Fig. 2: Two normalized histograms of the spatial distances between 
adjacent patches after the reordering, obtained for a noisy version 
of the image Barbara with unrestricted (a) and restricted (b) search 
areas. The lengths of the pathes in the spatial domain obtained with 
the unrestricted and restricted searches were 4.2 • 10^ and 6.1 • 10^ 
respectively. 

the overall computational complexity, and our experiments 
show that with a proper choice of B it also leads to improved 
results. The permutation applied by the matrix P is defined as 
the order in the found path. 

It is interesting to examine the characteristics of the patch 
ordering in the spatial domain. To this end we apply the patch 
ordering scheme described above to the patches of the noisy 
Barbara image shown in Fig. [Sj^a) with both unrestricted and 
restricted (B = 111) search neighborhoods, and with the 
parameters ^/n = 8 and e = 10^. We apply the obtained 
permutations to the patches, and calculate two normalized 
histograms of the spatial distances between adjacent patches, 
shown in Figure [2] Fig. |2ja) shows that only about 3% of 
neighboring patches in the path are also immediate spatial 
neighbors, and that far away patches are often assigned as 
neighbors in the reordering process. The histogram in Fig. 
|2jb) is limited to show only distances which are smaller or 
equal to 78, the maximal possible distance within the search 
window. It can be seen that despite the restriction to a smaller 
search neighborhood, only about 6% of neighboring patches 
in the path are also immediate spatial neighbors, and patches 
all over the search neighborhood are assigned as neighbors in 
the reordering process. 

In order to facilitate the cycle- spinning method mentioned 
above, we simply run the proposed ordering solver K times, 
and the randomness (both in the initialization and in assigning 
the neighbors) leads to different permutation results. We next 
describe how the quality of the produced images may be 
further improved using a subimage averaging scheme, which 
can be seen as another variation of cycle spinning. 

C. Subimage Averaging 

Let X be an n X (TVi - + l)(iV2 - + 1) matrix, 
containing column stacked versions of all the ^/n x ^/n 
patches inside the image Z. We extract these patches column 
by column, starting from the top left-most patch. When we 
calculated P as described in the previous section, we assumed 
that each patch is associated only with its middle pixel. 
Therefore P was designed to reorder the signal composed of 
the middle points in the patches, which reside in the middle 
row of X. However, we can alternatively choose to associate 
all the patches with a pixel located in a different position, e.g.. 



the top left pixel in each patch. This means that the matrix P 
can be used to reorder any one of the signals located in the 
rows of X. These signals are the column stacked versions of 
all the n subimages of size {Ni — ^Jn + 1) x {N2 — \/n + 1) 
contained in the image Z. We denote these subimages by Zj, 
j = 1, 2, . . . , n. An example for two of them, Zi and Z^, 
contained in a noisy version of the image Barbara, is shown 
Figlla). 

We already observed in fT2| and p3| that improved de- 
noising results are obtained when all the n subimages of a 
noisy image are employed in its denoising process. Here we 
use a similar scheme in order to improve the quality of the 
recovered image. In order to avoid cumbersome notations we 
first describe a scheme which utilizes a single ordering matrix 
P. Let Zj = RjZ be the column stacked version of Zj, where 
the matrix extracts the jth subimage from the image z. We 
first calculate the matrix P using the patches in X and apply 
it to each subimage Zj. Then we apply the operator H to 
each of the reordered subimages = Pzj, apply the inverse 
permutation P~^ on the result, and obtain the reconstructed 
subimages 

Yj = p-^HiPzj} = p-^i^{PR^z}. (6) 

We next reconstruct the image from all the reconstructed 
subimages fj by plugging each subimage into its original 
place in the image canvas and averaging the different values 
obtained for each pixel. More formally, we obtain the recon- 
structed image y as follows: 

n n 

y = ^ Rjy, = ^ Rjp-^i7{PR,z} (7) 

where the matrix Rj plugs the estimated jth subimage into 
its original place in the canvas, and 

n 

D = ^RjR, (8) 

is a diagonal weight matrix that simply averages the overlap- 
ping contributions per each pixel. When K random matrices 
P/e are employed, we obtain the final estimate by averaging 
the images obtained with the different permutations 

This formula reveals two important properties of our 
scheme: (i) the two summations that correspond to the two 
cycle- spinning versions lead to an averaging of nK candidate 
solutions, a fact that boosts the overall performance of the 
recovery algorithm; and (ii) if H is chosen as linear, then the 
overall processing is linear as well, provided that we disregard 
the highly non-linear dependency of P on z. 

D. Connection to BM3D 

The above processing scheme can be described a little 
differently. We start by calculating the permutation matrix P 
from the image patches x^. We then gather the patches by 
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Fig. 3: (a) Two subimages Zi and Z^ contained in a noisy version of 
the image Barbara, (b) Classification of the pixels in a noisy version 
of the image Barbara to centers of smooth (white) and non- smooth 
(black) patches. 



arranging them as the columns of a matrix in the order 
defined by P. This matrix contains in its rows the reordered 
subimages z^, therefore we next apply the operator H to its 
rows, and shuffle the columns of the resulting matrix according 
to the permutation defined by P~^. We obtain a matrix X, 
which contains in its rows the reconstructed subimages y^, 
and in its columns reconstructed versions of the image 
patches x^. We obtain the reconstructed image y from the 
patches Xj by plugging them into their original places in the 
image canvas, and averaging the different values obtained for 
each pixel. When K random matrices P/g are employed, we 
apply the aforementioned scheme with each of these matrices, 
and average the obtained images. 

Now, looking at the image processing scheme described 
above, we can see some similarities to the first stage of the 
BM3D algorithm. Both algorithms stack the image patches 
into groups, apply ID processing across the patches, return the 
patches into their place in the image, and average the results. 
We note that the BM3D algorithm also applies a 2D transform 
to the patches before performing ID processing across them. 
This feature can be easily added to our scheme if needed, and 
we regard this as a preprocessing part of the operator H. On 
the other hand, there are some key differences between the two 
schemes. First, while the BM3D algorithm constructs a group 
of neighbors for each patch, here we order all the patches 
to one chain, which defines local neighbors. Furthermore, 
this process is repeated K times, implying that our approach 
consider K different neighbors assignments. Also, while in the 
BM3D the patch order in each group is not restricted, ours is 
carefully determined as it plays a major role in our scheme. 
Finally, the ID processing applied by the BM3D consists of 
the use of a ID transform, followed by thresholding and the 
inverse transform, implying a specific denoising. Here we do 
not restrict ourselves to any specific ID processing scheme, 
and allow the operator H to be chosen based on the application 
at hand. We next demonstrate our proposed schemes for image 
denoising and inpainting. 

III. Applications and Results 

A. Image Denoising 

The problem of image denoising consists of the recovery of 
an image from its noisy version. In that case M = I and the 
corrupted image satisfies z = y + v. The patches x^ contain 



noise, and we choose the distance measure between x^ and Xj 
to be the squared Euclidean distance divided by n, i.e 

1„ 



(10) 



In our previous works fT2| and |T3| we applied a complex 
multi-scale processing on the ordered patches. Here we wish 
to employ a far simpler scheme; we choose a ID linear shift- 
invariant filter, and as we show next, we learn this filter from 
training images. Furthermore, we suggest to switch between 
two such filters, based on the patch content. 

We desire to treat smooth areas in the image differently than 
areas with edges or texture, as our experiments show that this 
approach leads to better results. More specifically, we employ 
different permutation matrices and filters in the smooth and 
non smooth areas of the image. We first divide the patches into 
two sets: Sg - which contains smooth patches, and - which 
contains patches with edges or texture. Let std(xi) denote the 
standard deviation of the patch x^ and let C be a scalar design 
parameter. Then we use the following classification rule: if 
std(x^) < C(j then x^ G Ss, otherwise x^ G Se- Fig [IJb) 
demonstrates the application of this classification rule to the 
noisy Barbara image shown in Fig. |3ja), where we use the 
parameters ^/n = 8 and C = 1.2 which we will later use 
in the denoising process of this image. White pixels are the 
centers of smooth patches and black pixels are the centers of 
patches containing texture or edges. It can be seen that the 
obtained image indeed contains a rough classification of the 
patches into smooth and non smooth sets. 

We next divide each subimage zj into two signals: zj^s - 
which contains the pixels corresponding to the smooth patches, 
and Zj e - which contains the pixels corresponding to the 
patches with edges and texture. We apply the nearest neighbors 
search method described above to the patches in the sets Ss 
and Se, and obtain two different permutation matrices P^ and 
Pg. Ps and Pg extract from zj the signals zj^s and zj^e, 
respectively, and then each apply a different permutation. We 
apply Pg and Pg to Zj^g and Zj^e and obtain z^^ and z^^, 
respectively, which are the signals to which we apply the 
filters. More formally. 













Z^ 




.^e{Zj,e} 







where we defined the matrix 
P = 



(11) 



(12) 



We next wish to find the filters hg and he applied to z*' ^ 
and z^g, respectively. We denote the convolution matrices 
corresponding to z^^, and z^^ by Z^^ and Z^^, and obtain 
the filtered subimages 



z^ 



7P 



-z> 



where we defined 



z^ 







z^ 



(13) 



(14) 
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The vector h stores the filter taps to be designed. We substitute 
( [T3] ) in ([7]), and obtain the reconstructed image 

n 

y = D-i^Rjp-iZ^h. (15) 

When K random matrices are employed, we obtain the 
final estimate by averaging the images obtained with the 
different matrices 



fe=i v j=i / 

= ^ED-i[Rr,...,Rj] 



/c=l 

where we defined 



h = Qh (16) 



Q = ^ED-MRf,...,Rj] 



k=l 



(17) 



Now let y^, j = 1, . . . , G be a training set which contains 
the column stack versions of G clean images. For each such 
image we create a noisy version by adding it noise with the 
same statistics as the noise in z. Then we calculate for each 
image z^ a matrix using ([17]), and learn the filters vector 
h by minimizing 

G 

h = argmin^ ||y^ - Q%f 

-1 



h 

G 



E(Q')V- (18) 



/c=l 



Once we have the filters vector h we can employ it to denoise 
z by building Q using ([17]) and then calculating 



y = Qh. 



(19) 



We can further improve our results by applying a second 
iteration of our proposed scheme, in which all the processing 
stages remain the same, but the permutation matrices are built 
using patches extracted from the first iteration clean result. 

In order to assess the performance of the proposed image 
denoising scheme we apply it to a test set containing noisy 
versions of the images Lena, Barbara and House, with noise 
standard deviations cr = 10, 25, 50. We learn the filters 
vector h from a training set containing the images Man, 
Peppers, Boat and Fingerprint. The parameters employed by 
the proposed denoising scheme for the three noise levels are 
shown in Table [U We note that the reason we chose a uniform 
filter length of 25 samples for all noise levels can be justified 
using Fig. [4] Fig. [4] shows the average of the PSNR values 
obtained in the first and second iterations for the 3 test images, 
as a function of the filter length, for the different noise levels. 
It can be seen that in both iterations the performance gain 
obtained using filters longer than 25 samples is negligible. 
The trained filters obtained in each iteration for the different 
noise levels are shown in Fig. [5] First, it can be seen that 



# Filter Length [samples] 



# Filter Length [samples] 



Fig. 4: Average of the PSNR values obtained for the 3 test images, 
as a function of the filter length, for the different noise levels: (a) 
First iteration, (b) Second Iteration. 





15 20 25 



Fig. 5: The trained filters learned from noisy images: Left column - 
the filters hs obtained in the first (top) and second (bottom) iterations. 
Right column - the filters he obtained in the first (top) and second 
(bottom) iterations. 



filters hs and hg indeed look different. It can also be seen 
that in the first iteration the shape of the filter does not 
change much as the noise level increases, and in the second 
iteration the filters obtained for the higher noise levels are 
similar, but very different from the filter obtained for a = 10. 
On the other hand, in both iterations the shape of the filter hg 
changes greatly as the noise level increases. 

For comparison, we also apply the K-SVD algorithm |[7|. 
The PSNR values of the results obtained with this algorithm, 
and two iterations of our denoising scheme are shown in Table 
[n] The noisy and recovered images obtained with our scheme 
for cr = 25 are shown in Fig. [6] First, it can be seen that the 
second iteration improves the results of our proposed scheme 
in all the cases but one. It can also be seen that the results 
obtained with two iterations of our scheme are comparable to 
the ones of the K-SVD for a = 10, but are much better than 
the ones of the K-SVD for a = 25 and a = 50. 

B. Image Inpainting 

The problem of image inpainting consists of the recovery 
of missing pixels in the given image. Here we handle the case 
where there is no additive noise, therefore v = 0, and M 
is a diagonal matrix of size N x N which contains ones 
and zeroes in its main diagonal corresponding to existing 
and missing pixels, correspondingly. Each patch may contain 




PSNR=32.48 dB 



PSNR=32.65 dB 



Fig. 6: Denoising results (PSNR) for the images Lena, Barbara and House (a = 25, input PSNR=20.18 dB): Left column - noisy images. 
Center column - 1 iteration results, and Right column - 2 iterations results. 
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Iteration 


K 




B 


C 


e 


Filters Length 


10 


1 


10 


6 


111 


1.6 


10^ 


25 


2 


10 


4 


441 


0.8 


10*^ 


25 


25 


1 


10 


8 


111 


1.2 


10^ 


25 


2 


10 


4 


441 


0.4 


10^ 


25 


50 


1 


10 


12 


111 


1.1 


10^ 


25 


2 


10 


5 


441 


0.2 


lO*^ 


25 



TABLE IL Denoising results (PSNR in dB) of noisy versions of 
the images Lena, Barbara and House,obtained with the K-SVD 
algorithm and the proposed scheme. For each image and noise level 
the best result is highlighted. 



Image 


Method 


cr/PSNR 


10/28.14 


25/20.18 


50/14.16 


Lena 


K-SVD 


35.49 


31.36 


27.82 


proposed (1 iter.) 


35.33 


31.58 


28.54 


proposed (2 iter.) 


35.41 


31.81 


29.00 


Barbara 


K-SVD 


34.41 


29.53 


25.4 


proposed (1 iter.) 


34.48 


30.46 


27.17 


proposed (2 iter.) 


34.46 


30.54 


27.45 


House 


K-SVD 


36.00 


32.12 


28.15 


proposed (1 iter.) 


35.83 


32.48 


29.37 


proposed (2 iter.) 


35.94 


32.65 


29.93 



missing pixels, and we denote by Si the set of indices of 
non-missing pixels in the patch x^. We choose the distance 
measure between patches and Xj to be the average of 
squared differences between existing pixels that share the same 



location in both patches, i.e. 



^(Xi,X^) 



{x,[k]-x,[k]y 



\s.ns,\ 



(20) 



We start by calculating the matrix P according to the 
scheme described in Section |n-B| with a minor difference: 
when a patch does not share pixels with any of the unvisited 
patches, the next patch in the path is chosen to be its nearest 
spatial neighbor. We next apply the obtained matrix to the 
subimages z\ and observe that the permuted vectors = Pz-^ 
contain missing values. We bear in mind that the target signals 
= Py^ should be smooth, and therefore apply on the 
subimages z^ an operator H which recovers the missing values 
using cubic interpolation. We apply the matrix P~^ on the 
resulting vectors and obtain the estimated subimages fj. The 
final estimate is obtained from these subimages using We 
improve our results by applying two additional iterations of 
a modified version of this inpainting scheme, where the only 
difference is that we rebuild P using reconstructed (and thus 
full patches). 

We demonstrate the performance of our proposed scheme 
on corrupted versions of the images Lena, Barbara and House, 
obtained by zeroing 80% percents of their pixels, which are 
selected at random. The parameters employed in each of the 
three iterations are shown in Table [nil In order to demon- 
strate the advantages of our method over simpler interpolation 
schemes we compare our results to the ones obtained by the 
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TABLE IIL Parameters used in the inpainting e xperiments 



Iteration 


K 


Vn 


B 


e 


1 


10 


16 


9 


10^ 


2 


10 


8 


43 


10^ 


3 


10 


5 


55 


10« 



neighborhood. Therefore, as calculating each of the distance 



TABLE IV: Inpainting results (PSNR in dB) of corrupted versions 
of the images Lena, Barbara and House with 80 percents of their 
pixels missing, obtained using triangle-based cubic interpolation 
(tri. cubic), overcomplete DCT dictionary, 1 and 3 iterations of the 
proposed scheme. For each image the best result is highlighted. 



Image 


Method 


PSNR 


Lena 


tri. cubic 


30.25 


DCT 


29.97 


proposed (1 iter.) 


30.25 


proposed (2 iter.) 


31.8 


proposed (3 iter.) 


31.96 


Barbara 


tri. cubic 


22.88 


DCT 


27.15 


proposed (1 iter.) 


27.56 


proposed (2 iter.) 


29.34 


proposed (3 iter.) 


29.71 


House 


tri. cubic 


29.21 


DCT 


29.69 


proposed (1 iter.) 


29.03 


proposed (2 iter.) 


32.1 


proposed (3 iter.) 


32.71 



measures (10) and (20) requires 0{n) operations, the number 
of operations required to calculate a single reordering matrix 
Pfe can be bounded by 0{NB'^n). Next, applying the matrices 
P/e and P^^ to the n subimages zj require 0(nN) operations, 
and so does applying either one of the operators H described 
above to the n subimages tF-. Finally , constructing an estimate 
image by averaging the pixel values obtained with the different 
subimages also requires 0{nN) operations, and averaging the 
estimates obtained with the different matrices P^ requires 
0{KN) operations. Therefore when K permutation matrices 
are employed, the total complexity is 

0{{n + K)N) + K [0{NB^n) + 0{nN)] = 0{NKB'^n) 

(21) 

operations, which means that, as might be expected, the overall 
complexity is dominated by the creation of the permutation 
matrices. For a typical case in our experiments, N = 512^, 
K = 10, n = 64 and B = 111, and the above amounts to 
1.86 • 10^^ operations. We note that while in our experiments 
we employed exact exhaustive search, approximate nearest 
neighbor algorithms may be used to alleviate the computa- 
tional burden. 



matlab function "griddata" which performs cubic interpolation 
of the missing pixels based on Delaunay triangulation p8| , 
|T9| . We also compare our results to the ones obtained using 
the algorithm described in chapter 15 of 1 16], which employs 
a patch-based sparse representation reconstruction algorithm 
with a DCT overcomplete dictionary to recover the image 
patches. We use a patch size of 16 x 16 pixels in order to 
improve the results this method produces. We note that we 
do not employ the K-SVD based algorithm which was also 
described in this chapter, as our experiments showed that it 
produces comparable or only slightly better results than the 
redundant DCT dictionary, at a higher computational cost. 
The PSNR values of the results obtained with the different 
algorithms are shown in Table IV Fig|7] shows the corrupted 
and the reconstructed images, with the corresponding PSNR 
values, obtained using triangle-based cubic interpolation, over- 
complete DCT dictionary, 1 and 3 iterations of the proposed 
scheme. First, it can be seen that the second and third iterations 
greatly improve the results of our proposed algorithm. It can 
also be seen that the results obtained with three iterations of 
our proposed scheme are much better than those obtained with 
the two other methods. 

C. Computational Complexity 

We next evaluate the computational complexity of a single 
iteration of the two image processing algorithms described 
above. We note that for the image denoising scheme, we 
assume that the filters training has been done beforehand, and 
exclude it from our calculations. First, building the matrix X 
which contains the image patches requires 0{nN) operations. 
We assume that when the nearest-neighbor search described 
above is used with a search window of size B x B, most of 
the patches do not require to calculate distances outside this 



IV. Conclusions 

We have proposed a new image processing scheme which is 
based on smooth ID ordering of the pixels in the given image. 
We have shown that using a carefully designed permutation 
matrices and simple and intuitive ID operations such as linear 
filtering and interpolation, the proposed scheme can be used 
for image denoising and inpainting, where it achieves high 
quality results. 

There are several research directions to extend this work 
that we are currently considering. The first is to make use 
of the distances between the patches not only to find the 
ordering matrices, but also in the reconstruction process of 
the subimages. These distances carry additional information 
which might improve the obtained results. A different direction 
is to develop new image processing algorithms which involve 
optimization problems in which the ID image reorderings act 
as regularizers. These may both improve the image denoising 
and inpainting results, and allow to tackle other applications 
such as image deblurring, where the operator M is no longer 
restricted to be point-wise local. Additionally, the proposed 
image denoising scheme may be improved by dividing the 
patches to more than two types, and treating each type differ- 
ently. Finally, we note that in our work we have not exhausted 
the potential of the proposed algorithms, and the choice of 
different parameters (e.g., B^e) for each set of patches may 
also improve the produced results. 
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